# Required packages ####
#install.packages("quantreg")
#install.packages("orthopolynom")
#install.packages("plotrix")
#install.packages("doParallel")
require(quantreg)
require(orthopolynom)
require(doParallel)

# Folders ####
directory<-Sys.getenv("US_Ineq_Repl")
CHwd<-file.path(directory, "Scripts/CH")
Datawd<-file.path(directory, "Processed")
Figwd<-file.path(directory, "Results/Figures")

# Functions #### 
setwd(CHwd)
source("CHfun.R")

beta1 <-function(t){
  b1<- 0.2 + 0.05*t
  return(b1)
}

beta0 <-function(t){
  b0<-0.5*rep(1,length(t)) 
  return(b0)
}

# Code begins ####
# Simulated data for regression ####
## Average age 35, with SD of 8 
m<-35
s<-8
cm<-cumsum(dnorm(20:60,mean=m,sd=s)/sum(dnorm(20:60,mean=m,sd=s)))
plot(20:60,cm)

N<-1000
set.seed(69)
un<-runif(N,min=0,max=1)
age<-matrix(rep(0,N),nrow=N,ncol=1)

for (i in 1:N){
  age[i]<-20+sum(cm<=un[i])
}

hist(age)
min(age)
mean(age)
max(age)

b0<- 0.5
b1<- 0.2

ly<-b0+b1*age
u <- runif(N)*(0.05*age) 
ly_sim=ly+u
plot(age,ly_sim)

min(ly_sim)
mean(ly_sim)
max(ly_sim)
age2=age^2

formula <- as.formula("ly_sim ~ age ")
## Simulated data
D<-as.data.frame(cbind(ly_sim,age,age2,rep(1,N)),ncol=4,nrow=N)
colnames(D)<-c('ly_sim','age','age2','earnwt')

## Save Simulated Data 
setwd(Datawd)
save(D,file="Data_Sim.Rda")

# Figure 1: Random Sample and Selected Conditional Quantile Functions ####
be01<-rq(formula,data=D,tau=0.01)$coef
be25<-rq(formula,data=D,tau=0.25)$coef
be50<-rq(formula,data=D,tau=0.50)$coef
be75<-rq(formula,data=D,tau=0.75)$coef
be99<-rq(formula,data=D,tau=0.99)$coef

age<-D$age
age2<-D$age2
ly_sim<-D$ly_sim

setwd(Figwd)
png('AppCData.png',width = 500, height = 500)
plot(age,ly_sim,type='p',pch=20,cex=0.6,col="dodgerblue2",xlab="X",ylab="Log Y")
lines(20:60,be01[1]+be01[2]*20:60,col=c("darkorange"), lty="solid",type="l",lwd = 2)
lines(20:60,be25[1]+be25[2]*20:60,col=c("darkorange1"), lty="solid",type="l",lwd = 2)
lines(20:60,be50[1]+be50[2]*20:60,col=c("darkorange2"), lty="solid",type="l",lwd = 2)
lines(20:60,be75[1]+be75[2]*20:60,col=c("darkorange3"), lty="solid",type="l",lwd = 2)
lines(20:60,be99[1]+be99[2]*20:60,col=c("darkorange4"), lty="solid",type="l",lwd = 2)

text(48.8, 0.3 +be99[1]+be99[2]*50,labels="99%",cex=1.3,col="darkorange4")
text(51.8, 0.1 +be75[1]+be75[2]*52,labels="75%",cex=1.2,col="darkorange3")
text(53.8,-0.4 +be50[1]+be50[2]*54,labels="50%",cex=1.2,col="darkorange2")
text(55.8,-0.5 +be25[1]+be25[2]*56,labels="25%",cex=1.2,col="darkorange1")
text(57.8,-0.7 +be01[1]+be01[2]*58,labels="1%",cex=1.2,col="darkorange")
dev.off()

# Estimates using various polynomial orders ####

setwd(Datawd)
load("Data_Sim.Rda")
formula <- as.formula("ly_sim ~ age ")
h<-c("(Intercept)","age")
repe<-1000
p<-length(h)
Q<-c(0,0.25,0.5,0.75,1)
q<-length(Q)
meth<-"pfn"
tot_cores<-detectCores()
no_cores <- max(5,round(0.25*tot_cores) )
no_cores <- 5*(no_cores%/%5)+if(no_cores%%5>2){5}else{0}
J<-69
taus <- 1:J/(J+1)

n<-2

CoefRq<-ParCoefRq(J=J,formula,Data=D,w="",mth=meth,nodes=1,ParRq=ParRq)

setwd(Figwd)
png('AppCint.png',width = 500, height = 500)
plot(taus,CoefRq[1,],type='p',pch=.20,cex=.6,col="dodgerblue2",xlab="Quantile",ylab="Bo")
lines(taus,beta0(taus),col=c("darkorange"), lty="solid",type="l",lwd = 2)
dev.off()

setwd(Figwd)
png('AppCslp.png',width = 500, height = 500)
plot(taus,CoefRq[2,],type='p',pch=.20,cex=.6,col="dodgerblue2",xlab="Quantile",ylab="B1")
lines(taus,beta1(taus),col=c("darkorange"), lty="solid",type="l",lwd = 2)
dev.off()

cls <- makeCluster(no_cores, type = "SOCK")
registerDoParallel(cls)

for (K in 2:10) {
  Mean<-matrix(rep(0,q*(p)),ncol=q,nrow=p)
  Variance<-matrix(rep(0,q*(p)),ncol=q,nrow=p)
  
  st <- Sys.time()
  results<-foreach(rr=1:repe,.combine='rbind',.inorder=F)%dopar%{
    Impact(n=K,nodes = 1,J=J,s=rr,formula = formula,h=h,Data=D,
           ParRq = ParRq,
           Legendre=Legendre,
           polynomial.values=polynomial.values)
  }
  et <- Sys.time()
  dis<-et - st
  cat("order = ", K, "\n")
  cat(paste0('time = ',round(dis,2),' ',dis%.%units),'\n')
  setwd(Datawd)
  volume<-paste0('_',repe,'_',K)
  name_file<-paste("sim_results_by_Q0",volume,".csv",sep="")
  write.csv(results, name_file)
  
  for (qq in 1:q){
    index<-seq(from = qq, to = repe*(q), by = q)
    Mean[,qq]<-colMeans(results[index,])
    Variance[,qq]<-colVars(results[index,])
    nam <- paste("Q", qq, sep = "")
    assign(nam,results[index,])
    setwd(Datawd)
    name_file=paste("sim_results_by_",nam,volume,".csv",sep="")
    fil=paste("write.csv(",nam,", name_file)")
    eval(parse(text=fil))
  }
  
  setwd(Datawd)
  name_file=paste("sim_results_mean_by_Q",volume,".csv",sep="")
  write.csv(Mean, name_file)
  
  name_file=paste("sim_results_variance_by_Q",volume,".csv",sep="")
  write.csv(Variance, name_file)
}

stopCluster(cls)




